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Introduction 

The fluid dynamics of flows of dispersed materials m a fluid is fundamental to susp< n 
sions, bubbly liquids, droplet flows, pneumatic transport, fluidization and erosion. Equa- 
tions of motion to describe these materials must deal with the interactions between them as 
well as the deformation of the carrier fluid. Models that treat assemblages of solid paitich s 
have been proposed and studied (Jenkins and Savage, 1983) that result in the particles 
behaving like a gas, with a pressure due to the fluctuations in the velocities that is at- 
tributed to colli sional motions of the indivual particles. Models for fluid-particle mixtures 
(Drew 1983) do not include this effect (Passman 1989) The purpose of the present paper 
is to derive constitutive equations to supplement the equations of motion that include the 
effects of the particle velocity fluctuations. The particle motions are assumed to be at a 
sufficiently high Reynolds number that the fluid motion is inviscid, but viscous effects such 

as boundary layer separation are neglected. 

Equations of Motion 

The appropriate general average is the ensemble average. The ensemble avc.ia.ge. is 
the appropriate generalization of adding the values of the variable for each realization, and 
dividing by the number of observations. We shall refer to a “process as the set of possible 
flows that could occur, given that the initial and boundary conditions are those appropriate 
to the physical situation that we wish to describe. We refer to a “realization” of the flow 
as a possible motion that could have happened. Generally, we expect an infinite number 
of realizations of the flow, consisting of variations of position, attitude, and velocities of 

the discrete units and the fluid between them. 

If / is some field (he., a function of position x and time t) for some particular real- 
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ization, //, of the process, then the average of / is 

(1) /(x,<)= / 

Jm 

where drn(fi) is the measure (probability) of observing process /i and M is the set of all 
processes. We refer to M as the ensemble. The ensemble average is the fundamental 
average that allows the interpretation of the phenomena in terms of the repeatability of 
experiments. Any one exact experiment or realization will not be repeatable; however, any 
lepetition of the experiment will lead to another member of the ensemble. 

In order to apply the averaging procedure to the equations of motion, we shall need 
some results about the averaging procedure. We shall also give these results for time- and 
v ol nine aver agi ng . 

In order to average to the exact equations, we need expressions for df/dt and V/. If 
f well behaved , then it is clear from the definition of the ensemble average that 

( 2 ) W = d l 

dt dt 

and 

(3) vj = V/ 

Functions are generally discontinuous at the interface in most multiphase flow. They 
are well behaved within each phase, however. Thus, consider WV/, where X k is the phase 
indicator function for phase k: 

(4 ) X k = / lf X G 

{ 0, otherwise. 

In the averaging process we will require the result 

( 5 ) + v, • VX k = 0. 

This relation has a reasonable physical explanation. Note that it is the “material” deriva- 
tive of X k following the interface. If we look at a point that is not on the interface, then 
either X k = 1 or X k = 0. In either case, the partial derivatives are both zero, and hence 
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the expression (5) is zero. If we consider a point on the interface, if we move with the 
interface, we see the function AT as a constant jump. Thus, its material derivative is zero. 

The averaged equations are 
Mass 

(C) + V • = J{ v - v, ) • VX k 


Momentum 

dX kpv 

(7) — ~~ + V • X k pvv = V • X k T + X k pg + (pv(v - v.) - T) • VX k . 


Next, we define the appropriate average variables describing multiphase mechanics. 
First, the geometry of the exact, or microscopic situation is defined in terms of the 
phase function X k . The average of AT is the average fraction of the occurrences of phase 
k at point x at time t. 

(S) a* = Xl 


The quantity dX^/dn is the delta function defining the interface, its average is the inter- 
facial area density. 


(9) 


ci 2 = 


dX k 

On 


All the remaining variables are defined in terms of weighted averages. The main, or 
“phasic” variables are either phasic weighted variables (weighted with the phase function 
AT) or mass- weighted (or Favre) averaged (weighted by X k p). 

The “conserved” variables are the density 


( 10 ) 


P k = X k p/a k , 


and the velocity 

(11) V x k p = X k pv/a k p k 
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The averaged stress is defined by 


( 12 ) 


T l = X k T /a k 


The average body force is 


(13) 


SV = X k pg/a k p x k 


As discussed above, several terms appear representing the actions of the convective 
and molecular fluxes at the interface. The convective flux terms are the mass generation 


rate 


(14) r* = p(v - Vi) ■ vat 
and the interfacial momentum flux 

(1 5 ) vsr* = pv(v - V,) • vx k 

The interfacial momentum source is defined by 

(1G) = -T • VAT 

The motion of the interfaces gives rise to velocities that are not “laminar” in general. 
The velocity fluctuations may be due to turbulence or to the motion in the phases due 
to the motion of the interfaces. The effect of these velocity fluctuations, whatever their 
source, on a variable is accounted for by introducing its fluctuating field (denoted by the 
prime superscript), which is the difference between the complete field and the appropriate 
mean field. For example, 

v * = v - V 

Then, 


(17) 


X k pvv =X k p(\ x k p + v' t )(v^ + v^) 

+ x^y k 
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The Reynolds stress is defined by 


(18) Tf e = -X k px' k v' k /a k , 

The averaged interfacial pressure p kl and shear stress r ki are introduced to separate 

mean field effects from local effects in the interfacial force. The interfacial piessme is 
defined by 

(19) Phi = pdXk/dnk/ai 

and the interfacial shear stress is 

(20) n-i = T k dX k /dn k /a t 

Thus, 

Mfc = — T • VX k 

=pVAT - r • VAT 

= Pkl XX k -T kl -XX~ k - T', t • V x k 

(21) =p kl Va k - r kt Va k + M' k , 

where we define the interfacial extra momentum source 

(22) M' t = M fc + p k iXa k - T kl ■ Xa k 


and introduce 


T' ki — —p' kl l + T ki — —(p — p k i) I + (t Tki)- 


The averaged equations governing each phase are 


Mass 


(23) 


dokPl 

dt 


+ V • a k p x k v x k 


p 


= T fc 


Momentum 


daknn p 

ot 


(24) 


+ V 


— x — xp — xp 

<*kp k v k p v k p 


= v • a k 


(t l + + a kPk g + M k + v£r* 
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The equation of conservation of mass for phase k (23) can be used in the momentum 
equation (24) to yield the Lagrangian form of the momentum equation: 


—X 

(*kp k 


- x D k vl P -r.(dvl P . 


Dt 


k — x 

= <*kPk 


(25) 


^-+*r-vvrj 

= V . a t (Tj + Tp ) + at -pls + Mi 
+ Pki Vat - Tkt ■ Va k + (v£ - v* p )T*. 


The jump conditions are derived by multiplying the exact jump condition by iij • VAj 
and averaging. This process yields the following conditions: 

Mass 


(26) 


rj + r 2 = o 


Momentum 


(27) Mj -f M 2 -f + v"'T 2 = m 

Constitutive Equations 

The exact equations of motion can be solved for the flow of an inviscid, incompressible 
irrotational fluid around an isolated sphere. We shall use this solution to derive information 
about constitutive equations for the force on the dispersed phase, the average stress, the 
Reynolds stress, and the interfacial pressure when the particle phase is allowed to have a 
random velocity. 

The fluid velocity at x for the irrotational flow of an incompressible inviscid fluid is 
expressed in terms of the velocity potential by 

(28) v ( x ) = V< i>(x) 

The continuity equation becomes 

(28) 0 = V • v = V • V<f> = V 2 <f>. 

The pressure at any point x is given in terms of the velocity by Bernoulli’s equation. 

(30) p - p{~^ + ^IW| 2 ) = Po = constant. 
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Consider a sphere located at a point z in a flow field, moving with velocity v ;) . The 
boundary condition at the surface of the sphere is 

(31) n • v p = n * v = n • V(f> at |x — z| = a, 

where a is the radius of the sphere, n is the normal to the surface of the sphere and v p is 
the velocity of the sphere. The boundary condition far from the sphere is 

(32) <f> -+ (jioo as |x - z| oo, 
where 

<f>oo = v 0 (f) * x + ■ e f • x 

is the velocity potential that would exist in the fluid if the sphere were not present. Here 
Vo(i) is the (unsteady) velocity of the fluid at the origin, and ef is the rate of strain tensor 
for the fluid. We shall assume that ej is constant. 

A convenient form for the solution of this problem is given by Voinov (1973), and is 

4> = v 0 (<) • x + • e/ • x 

+ ^ ( v p — v o(0 — z ' e /) • (x — z) 

(33) + l(x - z) • e 7 ■ (x - z) 

If there are many spheres in the flow field, the solution given above (33) will still be 
a good approximation if the spheres are sufficiently far apart that the flow disturbances 
due to the individual spheres do not interact. That is, the flow must be sufficiently dilute. 
Thus, we can think of each sphere as “isolated” in the sense that it only interacts with its 
neighbors through the averaged fields. We assume that each sphere lies in a “cell,” and 
inside that cell, the velocity is given by eq. (33). We shall approximate the cell to be a 
sphere of radius R. We choose R so that 

-7ra 3 /-7ri? 3 = a . 

3 / 3 
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Averaging 

W e now introduce a means of evaluating the averaging process. The first aspect of 
the ensemble average is that the sphere velocity is random, with a distribution function 
,/^( v ;»i x,f). The second aspect is that the sphere and the surrounding cell can lie with 
the sphere center anywhere within radius R of the space point x. The average is performed 
by first integrating over the distribution of the velocities that the sphere can have, followed 
by an integration over the possible positions that the sphere center can have. Note that if 
|x — z| < a, the material making up the sphere occupies the field point x and if |x — zj > a 
the fluid occupies the field point x. 

The average of a quantity z,v p ) is performed in two parts, that is first, we 

perform a conditional average of g for given sphere position z, integrating over the velocity 
space Vp, followed by the average over the spatial positions the sphere can have. We 
assume that the distribution of positions is such that 

dV j 

§7T R 3 _ X tVrf(x,<) 

is the probability of finding the sphere in a volume dV surrounding the point z, where 
x' = x - z. Thus, for the average over the fluid of a quantity g , we have 


flr(x,i|z)= / g(x, t; z, v„)/ (1) (v 7 j, z, t)dV v 


Here the notation </(x,f |z) is intended to suggest the conditional average assuming the 
sphere is located at z. The average of g over the fluid phase is then given by 


<7c(M) = 


r(f ? 3 - 


g(x, t\z)d£ldr, 


where 0(r ) is the sphere of radius r centered at x, and the integration is over the z variable. 

It will be convenient to introduce the average particle velocity and the fluctuation 
particle velocity as 



( 37 ) 


Vp(x,t) = Vp - Vp(x,<) 


no 



Note that 


(38) v' p = 0 

The particle kinetic energy per unit particle mass is 

(39) «?(*,*) = * / |Vp| 2 / (1) (v |M x,t) dV v 

~ Jn? 

and the Reynolds stress for the particles is defined by 

(40) T j e {x,t.) = -p d [ VpVp/ (1) (v p ,x,t) dV Vp 

Jr 3 

In order to evaluate the integrals appearing in the averaging process, we must express 
the z dependence of the velocities in terms of x and x f = x — z. We have 

v/(z) = v/(x) - x' • e/ 

and 

v p (z) = Vp(x) - x • e p 

where e p is the velocity gradient tensor for the average particle motion. We shall assume 
that this tensor is constant and symmetric. 

We have 


v(x, t; z, v p ) = V^(x, t; z, v p ) 

1 _ / a 3 

= v/(x) + “(v/(x) - Vp(x) - v' ( - x' • ( e/ - e p )) ( ^3 

- |( v /( x ) - Vp(x) - Vp - x' • (e/ - ep)) • x' f ^ 


(41) 


+i x ' e 4^r i x ' e/ x '(^) x ' 


Note that v/(x) is the fluid velocity that would exist at x if the sphere were not present, 
and vj(z) — \ p is the relative velocity between the sphere and the fluid evaluated at the 
sphere center. It is convenient to have expressions for the integrals of powers of x ! over 
S7(r). For these integrals, we note that 


(42a.) 


/ 

Jn(r) 


x' . . . x'dQ = 0 
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if the factor x' appears on odd number of times, and 


(42b) 


f d£l = 47r? -2 

JQ(r) 


(42c) 


/ 

JQ(r) 


x'xV/fl — -7T7 ,4 I 


(42d) 


/Q(r) 


XX XX (IQ = 7T7’ 6 S 

15 


where £ is a fourth order isotropic tensor defined in Cartesian coordinates by 


^ ijkl — “I” ^ik^jl "I - l&jk* 

We further note that if v is a vector, and e is a symmetric second order tensor with c tl = 0, 
then 

^ i j k / Vj C k l — 2 VjG j i . 

Derivation of Averaged Quantities 

In order to average eq. (41), we note that the average over the velocity fluctuations 
gives no contribution, by eq. (38). Then averaging over z gives 

(43) v* p (x,<) = v/(x,f). 


Similarly, to obtain the interfacial averaged velocity of the fluid, again the integration 
over the velocity fluctuations gives no contribution, and we have 


v c ,(x,f) 



Substituting and performing the integrations lead to the result 


(44) v ci (x,<) = v/(x,f). 

This result is a little surprising at first. The fluid at the surface of the sphere satisfies the 
condition n * v = n v p , but is allowed to slip in the tangential direction. After the passage 
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of the sphere, the fluid that was momentarily in contact with the surface of the sphere is 
again moving with the fluid. The result says that even during the time that it is m contact 
with the surface of the sphere, its average velocity is still equal to the average velocity of 

the fluid, and not of the sphere. 

Now let us compute averaged pressures using this formalism. The exact pressure can 
be computed by Bernoulli’s equation (30). In order to evaluate the derivatives m eq. (30), 
we note that x is constant during t derivatives, but dz /dt = v p = v p (z,t) + v p . Also, 
when evaluating V*, both t and z are held constant. The pressure is given by 


(45) 


(dv o 1 
p(x,t;z,v p ) = po - Pc l g t ' x + 2 


dv 0 dv p (x,t) 
dt dt 


+ (v p (x, t) + Vp) • e f - (Vp(x, t) + v'p) • e p 


■ x 



-i( v/ (x) - (Vp(z, t ) + vp - x' • (e/ - ep)) • (v p (z,<) + vp 
-|(v,(x) - (v,(z, t) + v' p )-x ■ (e, - e,)) ■ x' p) x' ■ (v,(i, t) + v p 


- ?(v„(z, 1) + vp • e, ■ x' (£) + |x' • e/ • x'x' . (v p (z, t) + V 1 ,) P ) + ,V/ • V/ 




+ ||v/(x) - (Vp(z,t) + vp - x' • (e/ e p)| y r6 
+ \v f • (v/(x) - (Vp(x, t) + vp - x' • ( e/ - e P )) ^- 3 
((v/(x) - (v,(x,f) + vp - x • (e/ - ep)) • x ) 

+ |(vp(x,<)+ v p-e/-x' (p) - ^(Vp(x,t) + vp-x / x'-e / -x' (p 


1 


(v/(x) - (v p (x,<) + Vp)) • e/ ■ x f 8 


[ (v/(x) - (v p (x, <) + vj,)) • x'x • e/ ■ x' 



a _ 8 

r io 


+ -(v/(x)-(vp(x,<) + Vp))-xx * e / ‘ x ( 10 
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wheic we have ignored terms of order e'j, e p , and eye p . Averaging over the velocity 
fluctuations gives 


p(x,t|z) = Po - Pc ( -x+ \ 

at 2 


dvp _ dvp^cj) 
dt dt 


+v ;) (x, t) - e f - v p (x, t ) ■ e p ] • x' ( - 


2 (v/(x) - v p (z, t) - x' • ( e/ - e p )) • v p (z, t ) - u* e (z, t ) 



-(vyCx) - v p (z ,t) - x' • (ey - e,,)) ■ x'x' • v p (z,i) + - — 


3x'.Tf(z,t)-x'l ( a 3 


A/ 



3 v /.(z, t) ■ e f ■ x' ( — ] + . e/ . x'x' ■ v p (z,t) f ^ ) + 


« 5 \ 1 


- , + _ V/ . V/ 


+ 


^l v /(x) - v p (z,t) - x' • (ey - e p )| + (z,f) 

+ 2 v / ‘ ( v /( x ) ~ v p ( x , 0 - x ' • fey - e p )) 





+ 


l ((v/(x) - Vp(x ,0 - x' . (ey - e ,j) • x ') 2 - 

8 Pd 


a 

T7o 


+ §v ;) (x, t) • ey • x' ( ~ ) - -v p (x,t) ■ x'x' .e r x'(l 



-fvy(x)- V p (x,t))-ey -X' 

-^(Vy(x) - V p ( X ,f)) • X'x' • ey • X' 

T^Mx) “ v p (x,<)) • x'x' • ey • x' 
The spatial integration is tedious, but results in 


( 46 ) 


—x <9v f 1 

Pc -Po- pc-^r ■ X - ^Vy(x) • vy(x) 


where we ignore terms of order a/R in addition to those ignored previously. We als 


ilso obtain 


( 47 ) 


p Ci (x, t) = p* - ~p c \y c (x,t) - v d (x,t)| 2 + \p c u l d u . 
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The interfacial force density M d is given by M d - p^^d- Thnt 


MX ' 


* |z) ( 1_ ^rr a 


W i* R* J Q( a) V <*d / 

This can be computed by substituting eq. (45) for the pressure, and recognizing that 

n = x'/a. We must also expand the terms uf e (z,t) = u ? e ( x >0 “ x * ( x ’^ <,iul 

Tf (z,t) = T «*(x,t) - x' • VTf (x,t) The result is that 

(49) Mj(x,t) = 

, -r(l f5S^_?S!^ + ^P( Xit ).VvJ^ 


+ «(iPc ( o 


2 dt 


_J[-(v^(x,f) - v7(x,t)) • (Vv^(x,t) - VV7(x,t)) J 

(?(v"(x,()-vj' , (x,i)) (v; , ’(x,( ) -v/lx.iliVaj 

+ l(v;»(x,() - v7(x,()(v;-(x,() - v7(x,f)) • Vai) 


7 ft e P c ^ mRe V7 rt j 7) X — u!} e — — —VOfrf’T^'' 
-a^-Vuf -a d — -V-Trf -* 2Q 

Note that no drag force is present ha eq. (49). This is the result of D’Alembert’s paradox, 
that is, there is no net force on a body moving at a constant velocity through an mvisad 

fluid at rest. 

If a distribution of stresses is applied to the surface of an elastic body, thcie results a 
distribution of stresses inside the body. These induced stresses are important in comput ing 
constitutive equations for solid-fluid mixtures. The average stress inside the sphere is given 


where 


t;(x,o=A/ // T < x -'w ,m 

av |7rcr Jo J Ju(r ) 

T(x,t|z) = [ T(x, t; z, Vp)/ ( 1 > (v, ( , z, t ) d\\ 

JR* 


Here T(x,t;z, v p ) is the stress at, point x inside a sphere having its center at z at time t. 
We shall assume that the spheres are linearly elastic solids, but we shall assume that the 
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deformation is sufficiently small that the fluid motion is unaffected by the deformation of 
the spheres. Then the stress-strain relation is given by 

( 5 °) T = g,[Vu + (Vu) <r ] + A S V • ul 

This can be written as 

( 51 ) T = <7 + 01 

where 

Ps'j V • u 

VulJ 

The spherical part of the stress satisfies (Love, 1932) 

(52a ) V 2 0(x, t; z, v p ) = 0 

®( x > z ! V P ) = ~p(x, t; z, v p ) on |x - z| = a 
Averaging over v p gives 

( 53a ) V 2 0(x,f|z) = 0 


= ^|[Vu + (Vu) tr ]-| 



(53b) 


0(x, t\z) = -p(x, t |z) on |x — z| = a 


Solving and performing the integration over z gives 

(54) e = 

The solution for <r is also given in Love (1932) and can be averaged in v p and then 
integrated in z to give 


°d(x,t) =p 


—2 

C 


20 


(v^(x, t) - V x /(x, t))(v*/(x, t) - V x /(x, t )) - i 


Re 


Pd 


(55) 


+ ^(|Vc P (x,f)-v7(x,t)| 2 +2^)l 
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We next turn to computations of the Reynolds stress using the velocity fluctuations 
due to the inviscid flow around a sphere. Using the expression for the velocity (41) and 
the average fluid velocity (43), we see that 


(56) 


v' c (x, t; z , v p ) = -(v/(z) - v p ) 



- |(v/(z) 

2 , 

H — x • er 
3 ; 




Averaging over the particle velocity fluctuations yields 


T? c (x,t|z) = -/0 c v' c (x,< |z)v' c (x,f|z) 


+ Pc 


1 ( a 


4 \ r 


*iRe 


3 (a 


Pd 


4 V r 


,8 


i Re 


|x'(x' . + (x' 

Pd 


-i /?e 

~ — IxM 




9 / a i i t / i 

+ 4 ) x x ( x ’ 


np Re 
• l 4 . 


Pd 


The integration over z can be performed, yielding 


Tf e (x,t) = ---o d p i 


20 


(W - v7)(vr - v7) 


rr Re 


(57) 


+3 ((v^ - v7) ■ (v- - v^) + 2 «J‘) I] - 


The fluid fluctuation kinetic energy is uf e = |v' • v' , and can be computed by taking 
the trace of eq. (57) for Tf e . The result is 


(58) 


= 7 a d 
4 


w - vn 2 + 


Re 


Conservation of Fluctuation Kinetic Energy 

In analogy with statistical mechanics for assemblages of molecules, the theory of av- 
eraging as applied to multiphase flows allows the computation of averaged equations for 

higher moments of velocity and pressure correlations. 

We start with the derivation of averaged equations for the fluctuation kinetic energy 
for each phase. The exact momentum equation is 

(59) + v • pvv = V • T + pg 
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We shall derive an equation for the evolution of the kinetic energy. The Lagrangian 
form of the momentum equation is 


(60) 


<9v 

P \ ^- + v • Vv ) = V • T -f pg 


If we take the dot product of v with eq. (60), we have 

1 


(61) 

We note that 


d\v 2 


dt 


+ v • V -v 2 1 = V ■ ( V • T) + px - g 


V • (V • T) = V • (T • v) - T : Vv 
If we also return to the Eulerian form, we have 


aw 


( 02 ) 


1 


dt + V • ^P v v = V • (T • v) - T : Vv + pv • g 
If we apply the ensemble average to eq. (62), we have 


dXkjpv 2 


1 


dt 


V • X k -pv 2 x = V • X*T • v - X k T : Vv 


(63) 


+ X : k p\ ■ g - [p-e 2 (v — v,-) + T • v) • XX k ] 


We define the fluctuation velocity of phase k by 


(64) 


Then 


v', = V - V- 


(65) 


yi — ( v k) 2 + 2v' fc • x c k p + {v x k p ) 2 


so that, noting that X^p\ f k = X kpv ~ X^p\ x k p — 0, we have 


A V = X„pK) 2 + X k p2v’ k ■ y’P + Xifaiy = a k ~p‘ k ut + a k -pl \ 


xp 


Furthermore, 



so that 


A^tPv =X k ph v' k )W k + X k p\(v' k fv x k p + v- • X k px' k V k + A>i(u/) 2 v/ 


jy . —x Re — ^P 

=Ot k q k + a kPk u k X k 


^•T 


Also, note that 

T v = T v/+ = T • v' t . 


Note further that 

T : Vv = T : Vv^ + T : Vv' fc 


Then we have 

A fc T • v = Ayr • v/ + X k T • V k = a k T x k ■ v/ - a k q‘ k 


and 


Ayr fVv = X k T : Vv/ + X k T : Vv^. = a k T k : Vv/ + D k , 


where q’f = q£ + q[. and D k = A*T : Vv' fe . 
Next, in the interfacial terms, we have 


>h 2 (v - Vi ) • va* = i(tf/) 2 r fc + v/ • v^r fc + p|«) 2 (v - v,o ■ VAy 


and 


(T • v) • VA'fc = (T • v/) ■ VA fc + (T • v'/ • VA t = M fc • v/ + Ky 

The equation for the conservation of fluctuation kinetic energy then becomes 


d -x 

dt akPk 


up + l(uP) 2 ) +V • c k -pl («P + /n 2 ) vP - V . ( ai v- . Tp) 


- V • a k {q k + qjf) + V • a*T fc • v/ + oc k p k v k ■ g 

- a k T X k : Vv/ - D* + + v/ • v^T* 

+ P-K) 2 (v - v , ) • VAfc + Mr- • v/ + W k 




(66) 
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The equation for the average kinetic energy is 


da kP% 2K P ) 
dt 


v* p • V • a k (T k + Tf e ) + M k ■ v xp 


(07) 


+ QkPk g ■ vr + v* p ■ v^Tfc. 


Subtracting this from eq. (66), we have 


da k pl 


:X Re 


dt 


( 68 ) 


— +V • a k plu^v x k P = ^ = Vv£' - V • a 4 (q* + qf) 
+ \(vl p ) 2 r k - ^K) 2 (v - V,) • VX k + W k - D k 


This equation has some interesting interpretations. First, note that the dissipation 
due to the Reynolds stress ct*Tf £ : Vv£ p acts as a source of fluctuation kinetic energy, 
while its counterpart for the molecular dissipation ctjtTj. : Vv^ does not appear in this 
equation. Dissipation on the macroscopic scale, then, winds up as different things on the 
microscopic scale. Also, the dissipation due to microscopic velocity fluctuations D k implies 
a loss of fluctuation kinetic energy. Thus, loss mechanisms, such as inelastic collisions or 
viscous dissipation in the velocity fulctuations, cause a loss of fluctuation kinetic energy 
to heat. Finally, the working of the fluctuations at the interface, W k appears as a source 
of fluctuation kinetic energy. 

Since this equation is unnecessary for the fluid phase, we shall ignore it for k = c. 
For k — d, we note that Dd = 0 is consistent with the linear elasticity assumption and the 
assumption that the particle radius a does not change. Furthermore, if assume no phase 
change (T* =0), and we ignore triple correlations in the particle velocity fluctuations 

( v p v p v p = °)> we have 

(69) 


du 

^ it +adrptd ' Vu ^d p = adT * e : + w * 


a dpd 


(70) 


Discussion of the Force on a Sphere 

The equations of motion for the mixture are 

da d 


dt 


+ V • a d v x / = 0 
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(71) 


da c 

Yk 


+ V • a c \ xp = 0 


(72) a d p x d ^ 




dt 


+ v7-Vv^ = 


— — (v c (x,t) — v d (x,t))(v c (x,t) — v d (x,t)) 
20 


-ojV(pJ + l?J«f - jSlv- - v7 1 2 ) + V ■ «.T 


d 


+v • < P 


+ 2 5 |v*' , (x,()-v;' , (x,()| 2 I 


+ a dPc ( 9 


^(M) _ + V^(x, t) • Vv* p (x, t) - v7(x, t) • Vv7(x, t) 


dt 


dt 


_1 [v7(x,0 - v7(x,0] • [Vv7(x,<) - Vv7(x, t)]j 
+p x c 0(v7(x, t) - v7(x, t)) - (v-(x, i) - v7(x, <))Vor d 
+ ^(v7(x, t) - v7(x, t))(v7(x, t) - v7(x, 0) • V«dj 


+/^V(« (i uf) 


(73) 


dt 


a cP*c l 


+ v7 ■ Vv7 ) = -a c Vp 


c 


1 


+V • < cc d p 

’ 20 


— X 

c 


(V-- - v;')(v>' - v 1 /) 


rpfie 

X r7 

Pd 


-oc d p c 


+3 ((v7 - v7) • (v7 - v7) + 2 uf) i] } 

+ (y*|v 7 (x,t) + v7(x,<)l 2 + Va 

0v7(x,f) 5 v7(x7) 


at 


at 


+ v7(x, t) • Vv7(x, t) - v7(x, t) ■ Vv7(x, t) 


_ — [v 1 ^ 
20 


[v‘0(x,i) - v7(x,()] ■ [Vv7(x,0 - Vv7(x,*)]) 
-p’ c Q(v7(x, t) - v7(x, ()) • (v7(x,<) - v’/(x,t))V aj 
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+ ^( v ?' , ( x . < ) - v 7( x - t ))( v ! f ( x < t )- v 7( x ’ t )) ■ Vo <i) 

+Pc‘^V( aj «"') + - («,Tf ) 

It is also possible to calculate the force on a sphere at z by computing 



where the integration is over the variable x', with x = z + x'. This results in 


(74) 


F p (z) 




dv 

dt 


f 1 

+ v / • e / + o 


dv f dv. 


dt 


dt 


+ vj ■ e f 


Note that this force agrees with Taylor’s (1928) calculation of the force necessary to hold 
a. sphere at rest in an accelerating stream, obtained by setting = 0 and v p = (). The 
force is 


(75) 


3 / 4 


— -7T a 


3 \ -x 


p x c Vf - Vv/ 


If we first take the gradients involved in eq. (72), using V-v^ p - ( 1 / )( da^/dt + v'/' 


Vtv ( /) and V ■ v 2 c p = — (1 /a c )(dad/ dt + v x c p -Vad) , then set v x d p 


0, u d c = 0, orf = const., 


and T d = 0; and assume one-dimensional, steady flow, then eq. (72) reduces to eq. (75) 
in the limit as a j —> 0. Moreover, it is clear that it should. Consider the one-dimensional 


situation pictured in Fig. 1. The continuum model for the particles between x and x + Ax 
gives the rate of change of the momentum of the particles and parts of particles between 
x and x + Ax, denoted by P^(x, .r + Ax), as the stress force transmitted to the particles 
by the particle parts outside of the interval, denoted by (o^i • TJ) \ x + W-i) TJlUa.. 
plus the force transmitted to the particles through their interface, denoted by M ,/ A x . 
Thus, 


(76) P d (x, x + Ax) = (a d i ■ T* d ) \ x + (a d (- i) • T^) | x+Ar + M d Ax 


The sum of the forces on all particles with their centers in the interval from x to 
./■ + Ax is equal to the sum of the pressure forces on all the particles involved. This is 
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denoted by E f pndS. This is equal to the rate of change of the momentum of all the 
particles, denoted by P p . Thus, 

(77) P p = E j pndS 

We note that eqs. (76) and (77) differ in the way they treat the particles being cut by the 
surfaces at x and x + Ax. The relation is that the 


P d(x,x + Ax) = P p + (a d i • TJ) U + (a rf (-i) * T x d ) | x +a* 


/ 

:ut,jn / 


pndS | x + £ cut, out / 


(78) 

Here 


^cut,in / pnrfS'lx-j-Ar d~ ^cut,out 


J pndSI 


i+ii 


:ut,in f pndS | j 


is the sum of the pressure forces on the surfaces of the cut particles at x whose centers are 
inside the interval from x to x + A.t, 


:ut,out J pndS\ j 
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is the sum of the pressure forces on the surfaces of the cut particles at x whose centers are 
outside the interval from x to x + Ax. A similar interpretation is valid for the cut particles 
at x + Ax. 

The terms on the right hand side of eq. (78) represent the resultants of forces on cut 
particles. If the approximate equation of motion inside the cut particle is V • T = 0, then 
the pressure force over the curved side, plus the stress resultant force over the flat side 
must add up to 0. (Note that if the particles are accelerating, then the forces add up to 
be the volume of the part of the cut particle, times the acceleration of its center of mass. 
Presumably, this force is small.) 

Conclusion 

Consistent forms for the interfacial force, the interfacial pressure, the Reynolds stresses 
and the particle stress have been derived for the inviscid, irrotational incompressible flow 
of fluid in a dilute suspension of spheres. The particles are assumed to have a velocity 
distribution, giving rise to an effective pressure and stress in the particle phase. The 
velocity fluctuations also contribute in the fluid Reynolds stress and in the (elastic) stress 
field inside the spheres. The relation of these constitutive equations to the force on an 
individual sphere is discussed. 
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